Effective thermostat induced by coarse-graining of SPC water 
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We investigate iiow tlie transport properties of a united atoms fluid with a dissipative particle dynamics ther- 
mostat depend on the functional form and magnitude of both the conservative and the stochastic interactions. We 
demonstrate how the thermostat strongly affects the hydrodynamics, especially diffusion, viscosity, and local 
escape times. As model system we use SPC water, from which projected trajectories are used to determine the 
effective interactions in the united atoms model. The simulation results support our argument that the thermostat 
should be viewed as an integral part of the coarse-grained dynamics, rather than a tool for approaching thermal 
equilibrium. As our main result we show that the united atoms model with the adjusted effective interactions 
approximately reproduce the diffusion constant and the viscosity of the underlying detailed SPC water model. 



I. INTRODUCTION 

In molecular dynamics, the addition of a thermostat is usu- 
ally motivated as a representation of interactions with the 
surroundings, important primarily for the thermal equilibra- 
tion of the system. The choice of thermostat generally de- 
pends on whether the system is characterized by constant 
temperature and volume (e.g. the Nose-Hooverii^ thermo- 
stat) or constant temperature and pressure (e.g. the Andersen^ 
or PaiTinello-Rahman'' thermostats). Under the standard as- 
sumption that the system is sufficiently chaotic (and therefore 
mixing), averages over the thermodynamic equilibrium en- 
semble can be calculated as time averages, provided that the 
system is allowed to equilibrate before measuring the time av- 
erage. This allows for measuring, e.g., the pressure and heat 
capacity of the system, and to investigate complex phenomena 
such as phase transitions. Due to the fluctuation dissipation 
theorem, thermodynamic properties defined through the par- 
tition function are unaffected by the choice of thermostat (for 
further discussion, see section Hlb. 

In contrast to the equilibrium properties, transport pro- 
cesses are intimately tied to the trajectories as they depend 
on the auto-correlation of velocities and forces through the 
Green-Kubo relations^. The thermostat changes how the 
trajectories approach equilibrium, and as a consequence the 
transport properties change as well. This is usually considered 
to be a problem; especially when running non-equilibrium 
(NEMD) simulations to measure transport properties, where 
the system is brought to a stationary, out of equilibrium, state 
by an external force field. The standard way of minimizing 
this effect is to make the system as large as possible, and the 
interaction with the thermostat as weak as possible. 

In this article we take a different view on the role of the 
thermostat: Rather than minimizing its effects, we view the 
thermostat as an integrated part of the dynamics. The ther- 
mostat is a mesoscopic representation of the exchange of en- 
ergy between the coarse-grained degrees of freedom and the 
microscopic, rapidly fluctuating, degrees of freedom. As a 
result we expect the effective coarse-grained dynamics to be 
of Langevin type. To demonstrate this point explicitly, we 
study a coarse-grained model of SPC water We use MD sim- 
ulations of SPC water as our microscopic system and define 



the projected dynamics as the center of mass motion (united 
atoms) of the SPC particles. We demonstrate that it is possi- 
ble to choose the effective conservative and stochastic inter- 
actions for the united atoms model, such that the equilibrium 
and transport properties are close to those of the projected mi- 
croscopic dynamics. 

The motivation for viewing the thermostat as a natural part 
of the system's dynamics lies in the Mori-Zwanzig theory on 
projection operators^i^iiSiiiii^. In short, the theory states that 
given a microscale dynamics, a lower dimensional represen- 
tation can be attained through a projection of the phase space 
(e.g. the map from the atomic coordinates to the center of 
mass of the molecules), where fast degrees of freedom can ei- 
ther give rise to Markovian (white) noise and dissipation, or 
be eliminated due to averagingii. Which of these two scenar- 
ios that best describe the system at hand depends on the ex- 
change of energy (heat) between the coarse-grained degrees 
of freedom and the degrees lost in the coarse-graining proce- 
dure. If the energy exchange is substantial, then the motion of 
the coarse-grained particles will not be smooth or determinis- 
tic, but rather evolve according to a stochastic (Langevin type) 
differential equation. We will show that our test system is de- 
scribed by the latter case. Naturally, how faithfully a particu- 
lar coarse-grained model represents the underlying dynamics 
depends strongly on the choice of projection, the mixing time 
of the fast degrees of freedom, and the details of the coupling 
between fast and slow dynamics. 

In principle, given a specific projection it is possible to 
derive the effective coarse-grained dynamics using the Mori- 
Zwanzig theory. In practice, however, the direct approach is 
neither computationally feasible, nor is the resulting coarse- 
grained system typically represented as a particle based sim- 
ulation. To make the method practically useful we use the 
united atoms coarse-graining as an ansatz. It is important that 
the forces in the coarse-grained model respects the fundamen- 
tal mechanical symmetries of the system, especially Galilean 
invariance, which implies conservation of linear and angular 
momentum. Therefore, we restrict ourselves to central forces 
that obey Newton's third law, and where the force between 
two particles depends only on the distance between the par- 
ticles. If we additionally assume that the stochastic compo- 
nent of the pair-wise forces is statistically independent, we 
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arrive at the dissipative particle dynamics (DPD) model, in- 
troduced in 1992 by Hoogerbrugge and Koelman— as a sim- 
ulation technique for complex hydrodynamic phenomena. In 
the DPD method, the thermostat is represented explicitly as 
dissipative and stochastic force a'^-''* . DPD can now be con- 
sidered a standard method for mesoscopic simulation. It has 
been used to study, for example, complex fluidsiS, sponta- 
neous self-assembly of amphiphilic molecules into bilayered 
membranesi^ vesicle s and hydrodynamics^^. 

We elaborate briefly on the direct consequences of choos- 
ing a Galilean invariant thermostat that faithfully characterize 
the local transport of linear and angular momentum. First, 
sufficiently close to equilibrium, the system obeys the classi- 
cal result of the asymptotic t^^l'^ decay of the velocity auto- 
correlation {d is the dimensionality of the system)^. In ad- 
dition, the interactions give rise to hydrodynamic modes in 
the fluid22i2i, which lead to the Navier-Stokes equations on 
the macroscopic level''*. Galilean invariance is especially im- 
portant for instance in NEMD measurements of the viscos- 
ity of complex fluid o'^i^^ . If the equations of motion are 
not Galilean invariant the thermostat causes a thermodynamic 
screening (see, e.g. Ref. [22 and references therein). A strik- 
ing example of the importance of correct hydrodynamic inter- 
actions is the block co-polymer melts studied by Groot and 
Madden^, which pass through a disordered meta-stable 
phase before organizing themselves in a hexagonal pattern. 
Using DPD allowed the system to go between the two regions, 
whereas a similar model without hydrodynamic interactions 
failed to pass through the barrier Hence, the local interactions 
between the coarse-grained molecules are important for the 
observed macroscopic state. Furthermore, the folding path- 
ways of many proteins and peptides are sensitive to the in- 
teractions with the surrounding water, and exhibit meta-stable 
intermediate states separated by kinetic barriers^^*^. 

In the standard approach to DPD the forces are usually not 
derived from the microscopic dynamics (for exceptions, see 
e.g. Ref. [23 and references therein). Instead, generic and sim- 
ple functional forms, such as constant or linearly decreasing 
up to a cut-off radius, are used. Moreover, the same func- 
tional form is usually used for both conservative and stochas- 
tic forces. It should be emphasized that this choice is guided 
by maximizing simplicity rather than strict physical argu- 
ments (the generic repulsive nature of the conservative force 
is, however, consistent with the effective interactions between 
the center of mass of clusters of particles^-). To account for 
transport properties it has become common practice to rescale 
time in order to match either diffusion or viscosity. A seri- 
ous and well documented pitfall of the DPD method in this 
version is that it cannot be tuned to match both the diffu- 
sion rate and the viscosity simultaneousl y' V^-^^i^' . In an at- 
tempt to overcome this problem, Junghans et al.— recently 
presented a study where a transversal component, orthogonal 
to the central force, was added to the DPD interactions. It 
was shown that this gives freedom to tune both the viscos- 
ity and the diffusion rate independently. However, it is clear 
that with this ansatz the conservation of angular momentum 
is no longer manifest, and we should therefore expect that the 
hydrodynamic behavior of the system is not necessarily faith- 



fully represented. In the current study we limit the interactions 
in the DPD dynamics to central forces but change the magni- 
tude and support of the interactions. We demonstrate that we 
can find parameter settings where both diffusion and viscosity 
are consistent between the coarse-grained and the microscopic 
model. 

The remainder of the article is organized as follows. First, 
we introduce DPD as a thermostat and show how the radial 
dependence of the stochastic force influences the approach 
to equilibrium. Second, we review how to estimate the con- 
servative force from the center of mass motion of SPC water 
molecules in atomistic MD simulations. Third, in the results 
section, we compare transport properties of the coarse-grained 
model to those of SPC water. Finally, we conclude with a dis- 
cussion and summarize the results. 



II. THE DPD THERMOSTAT 

In this section we review the DPD thermostat, and illus- 
trate how the structure of the stochastic force gives rise to 
Galilean invariance of the dissipative force. We also show 
how the structure of the noise affects the system's relaxation 
towards equilibrium. 

In its simplest form, the equations of motion for a DPD 
model, with particles positioned at r,, with velocities v,- and 
momenta p,-, can be written as a system of Langevin equations 
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F^ + Flj + F^ 



(1) 



where F^-, Fjj and F? are the conservative, dissipative, and 
stochastic forces between particles / and j. Both the conser- 
vative and non-conservative interactions in DPD are modeled 
by central forces obeying Newton's third law, ensuring that 
linear and angular momentum are conservedi^. In DPD, the 
stochastic force between particles / and j take the form 



Ff. = v/2fc^to(r,;)C,7e, 



(2) 



where r,^ is the distance between particles / and j, Cij is the 
unit vector pointing from j to ;, is Boltzmann's constant, 
and T is the temperature in Kelvin. The scalar function Co{rij) 
describe how the stochastic force depends on the distance be- 
tween the particles, and is interpreted as a symmetric Gaus- 
sian white noise term with mean zero and covariance 

{Qjit)Qf{t')) = {du>5jf + 5,f5j,)5{t-t'), (3) 

where 5ij and 5{t) are the Kronecker and Dirac delta func- 
tions, respectively. 

At thermal equilibrium the system is distributed according 
to the canonical ensemble 



/eq(x,p)=Z-'e-^('''P)ABr^ 



(4) 



where Z is the normalization term for the distribution. The 
equilibrium ensemble must be invariant under the equations 
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of motion. Since // is a constant of motion for Hamiltonian 
dynamics, the ensemble is invariant under this part of the dy- 
namics (this is true for any ensemble where the probability 
of finding the system in a given micro-state depends only on 
the energy), but in order for the dissipative and random forces 
to conserve the equilibrium, we need to choose the dissipa- 
tive force pp such that the combined contributions from the 
dissipative and stochastic forces cancel when acting on the 
equilibrium distribution: 



= J^Q-H{x,^)lk^T 



-FP + i^2A:BrA„(x)Vp^. 



-H{K.v,)lk^T 



-//(x,p)/iB7- 



(5) 



where ^ is the Fokker-Planck operator of Eq. ([T]), and A,y is 
a 3 X 3 matrix given by the covariance of the total forces on 
particles / and j. The equilibrium Fokker-Planck equation Q 
is commonly referred to as a fluctuation-dissipation relation. 
For the DPD model, this was first analyzed by Espanol and 
Warren—. Since it must hold for all points in phase-space, the 
only possible solution for the dissipative force is 



For tiie DPD model the force covariance is given by 



(6) 



- or (rij) e,/ ® e,7 when / 7^ j 
Y^k^i oP- {rik ) Cik <S) Cik when / = ; , 



(7) 



where Cg) denotes an outer product. Inserting Eq. (jTji into 
Eq. (|6]l we can write the dissipative force on a particle as a 
sum of pair-wise dissipative forces: 

Note that since depends only on the velocity differences 
between interacting particles, it is manifestly Galilean invari- 
ant, and it is clear from the derivation above that this is a di- 
rect consequence of the covariance property of the stochastic 
forces [c.f. Eq. (|3]l], which in turn stems from the assumption 
that Newton's third law applies. 

A subtle point is that there is no one-to-one relation be- 
tween the stochastic forces and A, the total force covariance 
matrix: Though F?- appears in the Langevin equation, the dy- 
namics depends only on A (as is seen from the Fokker-Plank 
equation for the system), and for each choice of the force co- 
variance matrix there are infinitely many choices of F? that 
obey Eq. (|7]i. For instance, the standard formulation of the 
DPD equations of motion contains N{N — l)/2 independent 
random variables, where is the number of particles, but it is 
possible to find a representation using no more than 3A^ inde- 
pendent random variables. In order to use this in the simula- 
tions, however, we would need to calculate the square root of 



the matrix A in each time step. Although both representations 
are equally valid, the standard form of the DPD equations of 
motion is much more efficient. 

The dissipation-fluctuation relation [c.f. Eq. Q] asserts 
that the systems thermal equilibrium is a fixed-point for the 
Fokker-Planck equation, or equivalently an invariant measure 
of Eq. ([1]). In order to better understand the effect of the 
stochastic forces on the path to thermal equilibrium it is illu- 
minating to study the time-evolution of the entropy of a non- 
equilibrium ensemble /(x,p) over the phase-space. Follow- 
ing Green—, we can express the difference in entropy of the 
ensemble / to the equilibrium distribution /eq as 



SH-S{t)= /dxdp/(x,p)log 



/(x,p) 

/eq(x,p) 



(9) 



which is negative for all ensemble distributions with the same 
support as the equilibrium distribution, and is zero if and only 
if the two distributions are equal. With strictly Hamiltonian 
dynamics, the entropy is constant in time. Intuitively, this is 
because the internal energy of the system needs to change in 
order for the ensemble to approach the equilibrium distribu- 
tion, but the Hamiltonian conserves the energy. With the addi- 
tion of dissipative and random forces it can be shown that en- 
tropy difference is a negative Lyapunov function on the space 
of ensemble distributions /: 

drS{t) ^k^T j dxdp/^ (^Vp, log A,j (^Vp^. log j-^ 
dxdp/^ 



e/i-(Vp,-VpJ log 



/ 



7. 



eq 



(10) 



which clearly is positive. In agreement with the second law 
of thermodynamics, the entropy continues to increase until 
f ~ feq- The main point of Eq. ( fTOl ) is that it shows ex- 
plicitly how the relaxation towards the equilibrium distribu- 
tion depends on the structure of the noise. As we have seen 
above, the conservative and dissipative forces follow from the 
Langevin equation and the shape of the equilibrium distribu- 
tion. Hence, for the purpose of capturing the transport prop- 
erties of the system, only the shape of the stochastic forces 
remains to be specified. 

Eqs. ([T]-[3] [8]) together establishes the general form of the 
DPD dynamics. Both the conservative force F^^, or equiva- 
lently the corresponding scalar potential, and the scalar func- 
tion (o{r) depend on the particular system of interest and need 
to be determined to obtain the correct DPD model. In prac- 
tice, this is the difficult part of DPD, and also the rationale 
behind the often used heuristic approach for deciding the in- 
teractions. In the following section we discuss a more system- 
atic approach for deriving the conservative interaction from a 
microscopic dynamics. 
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III. CONSERVATIVE POTENTIAL RECONSTRUCTION 

Several methods exist for deriving potentials from 
given equilibrium radial distribution functions (RDFs) 
27^34,35,36,37^38 (alternatively, direct time averagin g ov er the fast 
degrees of freedom can be used, see e.g. 



Refs. 



mi). The 



methods for reconstructing effective potentials from the RDF 
rely on the result by Henderson^", that two pairwise poten- 
tials resulting in the same RDF cannot differ by more than 
an additive constant. The importance of this theorem lies in 
the one-to-one correspondence between pairwise central force 
and the radial distribution function. 

In order to determine the effective potential corresponding 
to a given RDF, we use the inverse Monte Carlo method de- 
veloped by Lyubartsev and Laaksonen— . This method starts 
from a discretized Hamiltonian of the system. 



(11) 



which corresponds to using a stepwise constant potential, . 
Sa denotes the number of particle pairs separated by a distance 
in the range ra to ra+\, where ro—Q,ri ^dr and cx,- dr. 
The average of Sa is directly connected to the radial distribu- 
tion function, g{ra), by the relation 



(Sa) _ Vg 

N{N -l)/2 L^^ 



(ra), 



(12) 



where is the number of particles, L the side length of the 
simulation box and Va the volume of the spherical shell be- 
tween radii ra and ra+i. 



47Z 

Va = ^^ 

JO) 



'a+l 



(13) 



Using a start potential, , normally chosen as the poten- 
tial of mean force, 4>a* = -A:Bring(ra), a Monte Carlo (MC) 
simulation is made, where at each time step a random particle 
is moved a small distance in a random direction. From the 
simulation, the quantity (Sa) is measured, and the difference 
from the desired value, 5^ (given directly form the wanted 
RDF), is calculated. 



A{Sa) = {Sa) 



(14) 



Also measured is the correlation in the number of particle 
pairs at different distances in the system, {SaSy), and the co- 
variance matrix, {SaSy) — {Sa){Sy)- By taking these averages 
in the canonical ensemble, each entry in the covariance matrix 
can be related to the partial derivative of (Sa) with respect to 
the potential Oy, through 



d{Sa) _ {SgSy) - (5a) 



(15) 



As Lyubartsev and Laaksonen— suggests, this can be used in 
an expansion of {Sa), 

^(^«)=L3Sr^^*r + <5(A4>2), (16) 
7 ^^r 
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FIG. 1 : An outline of the procedure used for creating effective forces 
for a coarse-grained representation of SPC water. The radial distri- 
bution function is measured using the center of mass positions of the 
SPC molecules. With the inverse Monte Carlo technique, conserva- 
tive forces are then constructed, giving the correct RDF. By varying 
the parameters in the DPD thermostat, good agreement between the 
microscopic and coarse-grained dynamics can be obtained. 



from which an estimate of the error in the current potential, 

i.e. in is given by AOa- Updating the potential and 
restarting the MC simulation with the new potential will re- 
sult in a better approximation of the RDF, and after some iter- 
ations the method converges on a potential giving rise to the 
desired RDF. Difficulties with convergence in the method can 
normally be overcome by not changing the potential as much 
as described by AO, but instead chose the new potential as 
— - AA4>, where A € [0, 1]. For further discussion 
on the efficiency of the inverse Monte Carlo method, see e.g. 
Ref.lill 



IV. SIMULATIONS 



The procedure and technical details of our simulations are 
outlined below, corresponding to the schematics in Fig. [T] 
From a detailed MD simulation of SPC water, coarse grained 
conservative forces are constructed using the inverse Monte 
Carlo method. By also adding a DPD thermostat we wish to 
demonstrate how the thermostat can be tuned to obtain good 
agreement between the transport properties of the coarse- 
grained and the microscopic levels of description. It is im- 
portant to notice that the added thermostat does not affect the 
RDF, and therefore the conservative forces can be determined 
independently of the stochastic interactions. This is the ratio- 
nale behind decomposing the coarse-graining procedure into 
two steps, as shown in Fig.[T] 
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FIG. 2: Coarse-grained potential obtained by using the inverse MC 
method on center of mass data from an SPC simulation. Simulation 
results using 2180 particles, a temperature of 298 K, and a box with 
side length 4.06 nm. 



A. SPC water 

Simulations using the Simple Point Charge (SPC) model of 
water are performed using the molecular dynamics package 
Gromacs. The simulation is set up with 2180 water molecules, 
target temperature 298 K and target pressure 1 bar (NPT en- 
semble with Berendsen thermostat). Electrostatic interactions 
are modeled using the reaction field approach for separation 
distances larger than 1 .4 nm. 



B. United atoms 



As a coarse-grained description of the same system, we use 
spherically symmetric point particles representing the center 
of mass of the SPC molecules, and assume pairwise inter- 
actions between the coarse-grained particles. This type of 
coarse-graining is usually referred to as United Atoms (UA). 
The interaction potential between the particles is obtained by 
applying the inverse MC method ( section HIIli to RDF data for 
the center of mass of the SPC molecules. Figure |2] shows the 
coarse-grained potential as the piecewise constant approxima- 
tion given by the inverse MC method. In the coarse-grained 
simulations, the particle density and temperature are the same 
as in the original SPC simulation, but we now use the NVT 
ensemble (Berendsen thermostat), as this makes comparison 
with simulations using a DPD thermostat easier Simulations 
of SPC water with the NVT ensemble and box size 4.06 nm 
(which was the average box size using the NPT ensemble) do 
not show any differences compared to the NPT simulations, 
so we conclude that the choice of ensemble is not important 
for the analysis. 



C. United atoms with DPD-thermostat (UA-DPD) 

The UA-DPD simulations are performed using a standard 
velocity Verlet integration of the DPD equations of motion. 
This method gives second-order accuracy in the time-step for 
the conservative part, but is only first-order for the stochastic 
parti^. In the case of unbounded conservative interactions for 
small distances, the accuracy of the integration is mainly de- 
termined by the conservative part^; therefore, the integration 
error in the stochastic forces is of less concern. The simula- 
tion set-up is identical to the United Atoms simulation, with 
the addition of dissipative and stochastic forces. We examine 
the simple example where co{r) is linearly decreasing: 

fO, 0<r<ro, 
CO(r) = <^ (7(1 -r/rc), ro < r < r^, (17) 
[o, r,<r. 

The parameter a defines the strength of the dissipative 
force, and is the cutoff radius. Note that co{r) is also set 
to zero for r < ro- The reason for this is that for the rare occa- 
sions when the inter-particle distance goes below ro = 0.277 
nm, we do not want a stochastic force capable of pushing the 
particles closer together. This would result in uncontrolled 
magnitude of conservative forces, unless the time step is sig- 
nificantly lowered. 

In the simulations we vary the parameters a and to ex- 
amine the effect on the diffusion and viscosity values. Tuning 
the parameters, we try if it is possible to obtain diffusion and 
viscosity values consistent with those of the SPC simulation. 

D. Measurements of diffusion and viscosity 

The diffusion coefficients for SPC water and the UA model 
without DPD thermostat are obtained directly from Gromacs. 
In the UA-DPD simulations, the diffusion coefficients are cal- 
culated from the mean square displacement, using the Einstein 
relation 

Z, = limiM^^. (18) 

The viscosities for SPC and UA are measured with the NEMD 
method implemented in Gromacs^, while the viscosity of 
UA-DPD is measured applying a Poiseuille flow method^. 

V. RESULTS 

A. Self-diffusion and viscosity 

For the UA-DPD, we investigated the dependency of the 
diffusion and viscosity on the cutoff radius . By maintaining 
a constant a for different values of , we obtained the results 
shown in Figure [5] Both transport properties are strongly de- 
pendent on the value of r^: While the diffusion decreases with 
increasing r^, the viscosity increases. 
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FIG. 3: Diffusion and viscosity of UA-DPD simulations plotted as 
functions of the cutoff radius r^-. The strength of the dissipative force 
was identical for all simulations, a = 18.6 x 10^'^ kg m s^-'Z^. For 
the points where the error bar is not visible, the standard deviation is 
less than the width of the symbol. 



The diffusion of SPC water molecules was measured to 
D = 4.11 X 10^^ m^/s. Similarly, the viscosity was mea- 
sured to T] = 0.42 cP. As these values represents the cor- 
rect dynamics of the SPC water molecules and hence the 
coarse-grained representation, a coarse-grained model should 
be able to reproduce these numbers. The simplest approach, 
the UA model with a simple Berendsen thermostat, resulted 
in: D = 16.5 x 10"^ m^/s and rj =0.116 cP. Clearly, even 
though this model results in the same equilibrium distribution, 
i.e. the same RDF, as the SPC model, it does not represent the 
transport properties correctly. 

Adding the DPD thermostat, we tuned a separately for each 
Tf, such that the diffusion became approximately equal to the 
SPC value. The resulting cj values and the corresponding Kc 
values are listed in Table|I] Keeping D constant, we examined 
the variation of the corresponding viscosities as a function of 
rc and a. The results are summarized in Table HH The viscos- 
ity is at a minimum for small values of rc and then increases 
moderately with increasing r^.. The minimum value of the UA- 
DPD viscosity is close to the viscosity of the coarse-grained 
SPC system: 0.49 cP compared to 0.42 cP. In particular when 
comparing to the UA model, it is clear that the addition of 
a carefully tuned DPD thermostat plays an important role in 
preserving transport properties. 

We would like to stress that while our choice of weight 
function (o{r) was an appropriate starting point, it might be 
too simple to capture the system in the best possible way. Our 
results show that the transport properties depend on the func- 
tional form of (0{r). In the light of this it is expected, rather 
than surprising, that the viscosity is not matched exactly. In 
another study-^, we have suggested how the force covariance 
can be used to obtain an estimate of (o{r). A more elabo- 
rate theoretical framework based on this idea is currently in 
progress. 



TABLE I: The values of rc and a that results in a diffusion of ap- 
proximately D = 4.11 X 10^'' w?/s in UA-DPD simulations. This is 
under the condition that the simulation is set up according to Section 

rc [nm] G [IQ-" kg m s^^/^j 



0.36 


29.8 


0.38 


22.3 


0.40 


18.4 


0.50 


10.1 


0.60 


7.0 


0.70 


5.4 


0.80 


4.3 



TABLE 11: Diffusion and viscosity data for different simulations of 
water: SPC, UA and UA-DPD. The UA-DPD simulations are per- 
formed for different parameter values and CT of the linear weight 
function COi{r). Only the rc values are given in this table. The cor- 
responding values of CT are listed in Table |l] Data for real water is 
included as a comparison. All simulation measurements are done at 
298 K, the water data is valid for 298.2 K. 



Method 


rc [nm] 


D [10^^ m-s^i] 


T] [cP] 


SPC 




4.11 


0.42 ±0.004 


UA 




16.5 


0.1 16 ±0.007 


UA-DPD 


0.36 


4.16±0.07 


0.51 ±0.03 


UA-DPD 


0.38 


4.10±0.10 


0.49 ±0.02 


UA-DPD 


0.40 


4. 10 ±0.05 


0.49 ±0.02 


UA-DPD 


0.50 


4.07 ±0.07 


0.52 ±0.03 


UA-DPD 


0.60 


4.10±0.06 


0.52 ±0.03 


UA-DPD 


0.70 


4. 10 ±0.05 


0.55 ±0.03 


UA-DPD 


0.80 


4.08 ±0.07 


0.60 ±0.04 


Water 




2.27li5J 


0.8909li^l 



B. Escape time distribution 

In addition to determining the diffusion coefficient and the 
viscosity, we have also measured the escape probabilities for 
pairs of particles, that is, how large the probability is for two 
particles to separate beyond a given distance in a given time 
interval, given the initial separation. This gives a more de- 
tailed view of the local dynamics than the other transport prop- 
erties. The escape probabilities for SPC water, united atoms 
with the Berendsen thermostat, and the UA-DPD using a)(r) 
from Eq. ( fTTl ) are shown as level plots in Fig. ID The different 
levels in these figures represent the probability that two par- 
ticles, originally separated by the distance r, remain within a 
separation distance of 0.5 nm at a later time, f, displayed on 
the y-axis. 

Generally, the escape time is decreasing with increasing ini- 
tial separation r. Note especially that the average escape time 
is much larger in the region r < 0.3 nm than for r > 0.3 nm. 
This is consistent with the well-known caging effect. Parti- 
cles repeatedly bounce against their closest neighbors in the 
fluid, and for a pair of nearby particles to separate beyond 
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0.25 0.3 0.35 0.4 0.45 0.5 

r [nm] 



FIG. 4: Probability for a pair of particles to remain within a distance of 0.5 nm from each other during a time interval [0,t], given the distance 
r between the particles at time 0. The different figures show the escape time for (A) SPC water, (B) UA simulation with conservative force 
only, and (C) UA-DPD using 0}{r) from Eq. ^Ji with r,. = 0.38 nm and a = 22.3 x 10"'^ kg m s"^/^. 



the typical neighbor distance (about 0.28 nm for SPC water) 
requires either a concerted motion of several particles, or a 
very high kinetic energy of the particles. Note that this bar- 
rier is largely missing from the united atoms system with the 
Berendsen thermostat; this is the main reason why the diffu- 
sion rate is too high compared to the SPC water In contrast, 
the UA-DPD thermostat fitted to the diffusion rate and vis- 
cosity (Fig. HJZ) exhibit escape times similar to those of the 
SPC model. Hence, matching diffusion and viscosity seems 
to confer agreement with respect to more general transport 
properties. 



VI. SUMMARY AND DISCUSSION 

We have applied a coarse-graining scheme to the SPC wa- 
ter model. The effective dynamics of the center of mass of the 
water molecules is Langevin-like with a conservative (drift) 
component. The conservative part of the dynamics is repre- 
sented by a pairwise potential that can be determined from the 
radial distribution function. Stochastic and dissipative forces, 
motivated by the Mori-Zwanzig projection operator formal- 



ism, are also included as a DPD thermostat. We have demon- 
strated that the radial dependence of the stochastic and dissi- 
pative forces strongly affect the transport properties. The main 
conclusion is that the diffusion rate and the viscosity of the 
original system can be approximately matched by the coarse- 
grained model. This was obtained by tuning the magnitude 
and functional form of the random and dissipative forces. 

We argue that in meso- and microscale simulations of 
coarse-grained Hamiltonian systems, the stochastic and dis- 
sipative components of the interactions should normally be 
considered as an integrated part of the effective dynamics. 
For coarse-grained models with only pairwise interactions, a 
dissipative particle dynamics ansatz is well suited for repre- 
senting the non-conservative part of the dynamics, since it is 
the most general pairwise interaction that respects the local 
conservation or transport of linear and angular momentum, 
i.e. hydrodynamic behavior. The DPD interactions depend on 
the inter-particle distance. We argue that this functional de- 
pendence should the tuned so that the transport properties are 
consistent with the microscopic dynamics. More generally we 
show explicitly in Eq.[TO]how the relaxation towards equilib- 
rium depends on the structure of the thermostat. Finally, we 
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point out that all practical molecular dynamics methods, not 
only united atoms models, are ultimately coarse-grained from 
some more fundamental level, either from quantum mechan- 
ics or by removing fast degrees of freedom (e.g. vibration 
modes in covalent bonds) from a classical mechanics model. 
The interaction with these fast degrees of freedom plays an 
active role when the system equilibrate. The approach used in 
this paper can therefore possibly also be applied to atomistic 



molecular dynamics simulations. 
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